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Abstract 

We implement a new and accurate numerical entropic scheme to investigate the 
first-order transition features of the triangular Ising model with nearest-neighbor 
{Jnn) and next-nearest-neighbor {Jnnn) antiferromagnetic interactions in ratio R = 
Jnn/Jnnn = 1- Important aspects of the existing theories of first-order transitions 
are briefly reviewed, tested on this model, and compared with previous work on the 
Potts model. Using lattices with hnear sizes L = 30, 40, . . . , 100, 120, 140, 160, 200, 240, 360 
and 480 we estimate the thermal characteristics of the present weak first-order tran- 
sition. Our results improve the original estimates of Rastelli et al. and verify all the 
generally accepted predictions of the finite-size scaling theory of first-order transi- 
tions, including transition point shifts, thermal, and magnetic anomalies. However, 
two of our findings are not compatible with current phenomenological expectations. 
The behavior of transition points, derived from the number-of-phases parameter, 
is not in accordance with the theoretically conjectured exponentially small shift 
behavior and the well-known double Gaussian approximation does not correctly 
describe higher correction terms of the energy cumulants. It is argued that this dis- 
crepancy has its origin in the commonly neglected contributions from domain wall 
corrections. 
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1 Introduction 



Generalizations of the Ising model including further than nearest-neighbor in- 
teractions may give rise to complicated spatial orderings and produce complex 
and rich critical behavior [1,2]. The transitions between ordered and disordered 
phases may be continuous or of first-order with a tricritical point between 
them. The exactly soluble Ising model in d = 2 [3] with the addition of next- 
nearest-neighbor interactions becomes an analytically intractable problem and 
an understanding of the effects of adding such further couplings on the criti- 
cal behavior of the system is an open important problem of great theoretical 
and experimental interest. We will be particularly interested is cases involving 
competing interactions with ground-state arrangements that mimic a sublat- 
tice order or superantiferromagnetic (SAF) order in which ferromagnetic lines 
along the directions of the lattice alternate with lines of opposite oriented 
spins. Such models are of great theoretical and experimental interest. They 
simulate various types of antiferromagnets [4,5,6,7,8], but also important mod- 
els of alloy orderings [4]. Due to the well-known correspondence between spin 
systems and lattice gas, they provide an understanding of the behavior of ad- 
sorbed gases on crystal surfaces [9,10]. It should be noted here that, remarkable 
efforts have been made to predict the order of the transition according to the 
Landau-Lifshitz rules, which, on the basis of group theory arguments, select 
the ordered configurations that can be reached from the disordered phase by 
a continuous transition [4,11]. K. Binder [4] has reviewed these rules and has 
also pointed out well-known counter-examples, emphasizing the experimental 
difficulties in distinguishing a weak first-order transition from a second-order 
transition. 

Second-order transitions are more special than first-order transitions, but they 
are theoretically much better understood. At a first-order transition there is 
no diverging correlation length, and in general one cannot restrict attention 
to long wavelength phenomena, thus no universality, as in critical phenom- 
ena, is to be expected. Furthermore, it is well known that antiferromagnetic 
arrangements, with several order-parameter components, may be produced 
by the competition of the interactions [4]. However, no general theoretical 
agreement exists connecting the symmetry of spin structures, the number of 
order-parameter components, and the range of interaction with the expected 
critical behavior and in particular the order of the phase transition. Further- 
more, many authors have demonstrated the difficulties in properly identifying 
the order of the transition on a ground of high-temperature expansion, scal- 
ing, renormalization group transformations, and Monte Carlo simulations (see 
Ref. [4] and references therein). Roomany and Wyld [12] pointed out that the 
occurrence of second-order type effects at the weakly ffrst-order transitions 
can be explained by a comparison of the correlation length ^ with the lattice 
size L. In the case of the q — 5 Potts model's weak first-order transition. 
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Peczak and Landau [13] have observed pseudocritical exponents close to the 
conjectured values of the critical indices in the q — 4 Potts model. 

This paper is concerned with the analysis of numerical data, obtained via 
an accurate entropic Monte Carlo scheme, on triangular Ising finite systems 
with nearest-neighbor (Jn„) and ncxt-ncarcst-neighbor (Jnnn) antiferromag- 
nctic interactions in ratio R = Jnn/Jnnn = 1- Rastelli et al. [14] have recently 
presented numerical evidence (for R = 0.1, 0.5, and 1.0) that, in the ther- 
modynamic limit, this model undergoes a first-order phase transition, from a 
layered ordered phase (SAP phase) to a high temperature paramagnetic phase. 
Here, we will present a more detailed identification of the nature of this transi- 
tion by looking again at the size dependence of the traditional thermodynamic 
quantities but also by implementing the Lee-Kostelirlitz method [15,16] and 
detecting the first-order character of the transition via the size-dependence of 
the free energy barrier. Purthermore, we will improve the original estimates 
of Rastelli et al. [14] and try to verify the predictions of the finite-size scal- 
ing (PSS) theory of first-order transitions, including transition point shifts, 
thermal, and magnetic anomalies. 

The theory of PSS near second-order transitions has a rich and longstanding 
literature [17], starting with the pioneering works of Pisher [18] and Pishcr and 
Barber [19]. This theory has been extended to first-order transitions [4,15,16,20,21,22,23] 
and PSS and renormalization group methods have proven to be very useful 
tools for the study of first- and second-order phase transitions [15,16,17,18,19,20,21,22,23,24,25]. 
Our attempt here to elucidate the distinctive first-order features of the present 
model, by an extensive numerical study, will closely follow previous analogous 
studies carried out on the q = 5, q = 8, and g = 10 Potts model [15,16,21,26,27]. 
In these studies the existing theories of first-order transitions have been tested 
and verified but several important aspects have not been thoroughly clarified, 
especially in the cases of weak first-order transitions. Not only the demonstra- 
tion of a weak first-order transition is more difficult but also strong finite-size 
effects may obscure the true asymptotic behavior [16] and some theoretical 
predictions may not be met even at quite large lattices (since the correlation 
length may be very large). In such cases, it is very difficult to discriminate 
between wrong phenomenological expectations and correct theoretical predic- 
tions. Therefore, comparisons of different models and alternative studies may 
provide useful information and a better understanding of finite-size effects, 
something that is necessary for the correct interpretation of the numerical 
data and the verification of the theories. The present model, showing a weak 
first-order phase transition, offers the opportunity of such a contrasting test 
with the well-studied cases of the Potts model. 

Our interest in the present model stems from our recent studies on the anal- 
ogous square SAP model [28,29], which is again an Ising model on the square 
lattice with nearest- neighbor (J^n) and next- nearest- neighbor {Jnnn) antifer- 
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romagnetic interactions. The ground state of the square SAF model is four- 
fold degenerate and consists of the arrangements in which ferromagnetic rows 
(columns) alternate with opposite oriented spins. In the square model the SAF 
order can be obtained in both cases of a ferromagnetic or an antiferromagnetic 
nearest-neighbor couphng (for its T = phase diagram see Refs. [7,30,31]). 
Similarly, the present triangular SAF model, with antiferromagnetic interac- 
tions in ratio R = Jnn/Jnnn = 1, has a six-fold degenerate ground state and 
consists of the six arrangements in which ferromagnetic lines alternate with 
opposite oriented spins in the three lattice directions {T — phase diagrams 
are given in Refs. [5,9]). The Hamiltonian that governs these systems, in zero- 
field, is 

Jnn ^ ] (^i^^j ~l~ Jnnn ^ ^ ^i^ji (1) 

where here both nearest-neighbor (Jnn) and next- nearest- neighbor {Jnnn) in- 
teractions will be assumed to be positive (antiferromagnetic). This Hamilto- 
nian has been studied also in = 3 FCC lattices and evidence for both first- 
and second-order phase transitions have been presented [4,6], although the 
distinction of the order of the transition was difficult in that case too. 

The square model, governed by Eq. (1), develops at low temperatures SAF or- 
der for R > 0.5 and several early studies [7,30,31,32,33,34,35], mainly based on 
importance sampling [36,37,38,39,40], have suggested that this system possess 
anomalous exponents and a non-universal critical behavior with exponents de- 
pending on the coupling ratio it!. However, despite this early accepted scenario, 
there has been recently renewed interest and some attempts to re-examine the 
behavior of this model appeared. In several papers Lopez et al. [41,42,43] have 
used the cluster variation method (CVM) to study this model, concluding 
that the system undergoes a first-order transition for a particular range of the 
couphng ratio R {R — 0.5 — 1.2). On the other hand, this different scenario 
- predicting first-order transitions between ordered and disordered phases fol- 
lowed by continuous transitions outside the first-order region - has now been 
seriously questioned by the present authors, who have presented very strong 
evidence [28] for the case R — 1 that points out a well-obeyed second-order 
scahng behavior. The Buzano and Pretti [44] attempt to verify the scenario of 
Lopez et al. [41,42,43], by including an additional 4-body coupling, and using 
again the CVM, has revealed a further peculiar inadequacy of the CVM. The 
limiting case (J„„ = 0), where the exact solution of the Baxter model [45] 
applies, was also predicted by the CVM to undergo a first-order transition, in 
contradiction to the exact solution. 

The above introduction describes a rather controversial situation on the class 
of models with antiferromagnetic ground-state arrangements. Therefore, it is 
of interest to follow closely the previous FSS analysis applied to the Potts mod- 
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els and study in more detail the triangular SAF model. In this case, we will 
insist that this system undergoes a genuine, but weak, first-order transition, 
as originally predicted by Rastelli et al. [14]. Having a six-fold degenerate 
ground state arrangement, this model will be shown to produce first-order 
characteristics that he between the q — b and g = 6 Potts model. The rest of 
the paper is organized as follows. In the next Section we outline an extensive 
entropic sampling program. This program goes beyond our previous prac- 
tice in other apphcations [46,47,48,49] and is based on (i) the Wang-Landau 
(WL) method [50], (h) our dominant energy restriction (CrMES) scheme [46], 
and (iii) a second-stage improvement that combines the WL [50], the Lee en- 
tropic [51,52], and the broad histogram (BH) [53] or transition matrix [54] 
methods. In Section 3 we shall review all necessary theoretical aspects of first- 
order transitions that are then used for the analysis of our numerical data. The 
free energy barrier method of Lee and Kosterlitz [15,16], used in the literature 
for the identification of a first-order transition, is discussed in subsection 3.1.1, 
while the size-dependence of thermal properties derived from the well-known 
double Gaussian approximation is summarized in subsection 3.2. Finally, some 
new transition points derived from the number-of-phases parameter [26,27] are 
presented in subsection 3.3. Section 4 presents the analysis of our numerical 
data. In subsection 4.1 we make use of several alternative methods for the esti- 
mation of the transition temperature presenting a detailed analysis of the new 
transition points derived from the number-of-phases parameter and a discus- 
sion on their conjectured exponentially small shift behavior. In subsection 4.2 
we explore our data for the energy cumulants and the magnetic susceptibility, 
emphasizing on the examination of the higher-order corrections predicted by 
the double Gaussian approximation. The characteristic values of the energy 
cumulants, i.e. the coefficient of the dominant diverging term of the specific 
heat and the limiting values of the second- and fourth-order energy cumulants 
are determined and the corresponding theoretical predictions are critically dis- 
cussed. This serves as a useful self-consistency test of our numerical data but 
also of the theoretical predictions. Finally, our conclusions are summarized in 
Section 5. 



2 Numerical approach 



Computer simulations based on Monte Carlo sampling methods have increased 
dramatically our understanding of the behavior of systems of classical statisti- 
cal mechanics. The Metropolis method and its variants were for many years the 
main tools in condensed matter physics and critical phenomena [36,37,38,39,40]. 
However, for complex systems, effective potentials may have competing min- 
ima, or a rugged landscape, becoming more pronounced with increasing system 
size. In such cases, the traditional methods become completely inefficient, since 
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they cannot overcome such barriers and cross from one basin to another in 
the state space. On the other hand, a large number of "generahzed ensemble" 

methods have been proposed to overcome such troubles [39,40,50,51,52,53,54,55,56,57,58,59,60,61,62,63] 
One important class of these methods emphasizes the idea of directly samphng 
the energy density of states (DOS) and may be called entropic sampling meth- 
ods [39]. In an enropic sampling method, instead of sampling microstates with 
probability proportional to e~^^ , we sample microstates with probability pro- 
portional to [G{E)]~^, where G{E) is the DOS, thus producing a "flat energy 
histogram". The prerequisite for the implementation of the method, is the 
DOS information of the system, a problem that can now be handled in many 
adequate ways. Over the last two decades, there have been a number of inter- 
esting approaches addressing this problem. The most remarkable examples are 
the Lee entropic [51,52], the multicanonical [55,56], the BH [53], the transition 
matrix [54], and the WL [50] methods. In particular, the WL algorithm [50] 
has been one of the most refreshing improvements in Monte Carlo simula- 
tion schemes and has been already applied to a broad spectrum of interesting 
problems in statistical mechanics and biophysics [64]. 

Recently, the present authors have introduced [46,47] a dominant energy sub- 
space implementation of the above entropic methods, called critical minimum 
energy subspace (CrMES) method. The main idea is a systematic restriction 
of the energy space with increasing lattice size. The (WL) random walk takes 
place in the appropriately restricted energy space and this restriction pro- 
duces an immense speed up without introducing observable errors. For the 
temperature range of interest, that is the range around a critical or, for the 
present model, a first-order transition point, this scheme may be used to deter- 
mine all finite-size thermal anomalies from the final accurate DOS of the WL 
scheme. A further simplification, for the determination of the magnetic finite- 
size anomalies, was suggested [47] by using the WL approach as a one-run 
entropic strategy. In this CrMES WL-entropic samphng method the magnetic 
properties are obtained by recording appropriate histograms in the high levels 
(or final levels) of the WL process. It was argued and numerically verified that 
these high levels, in which the detailed balanced condition is quite well satis- 
fied, give good approximations for the microcanonical estimators, necessary for 
the evaluation of other properties, besides the thermal ones, of the statistical 
system. The method was efficiently combined with the N-fold way [37,47,65,66] 
in order to improve statistical reliability, but also to produce BH estimators, 
from the same high levels, for an additional calculation of the DOS. 

Let us now outline the main ingredients and equations used in the above de- 
scribed implementation. The approximation of canonical averages, in a tem- 
perature range of interest, is as follows: 

^ EE{Q)EG{E)e-^^ ^ j:EeiE,,E,){Q)E,WLGiE)e-^^ 

j:EG{E)e-^^ EE,iE,,E.)G{E)e-^^ ' ^ ^ 
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The restricted energy subspacc (Ei. E2) is carefully chosen to cover the tem- 
perature range of interest without introducing observable errors. The micro- 
canonical averages (Q) e are determined from the {E, Q) -histograms (denoted 
by Hivl(E, Q)), which are obtained during the high levels of the WL process 

{Q)e = {Q)e,wl ^ E Q ^r^ff ^ Hwl{E) = E H^l{E, Q) (3) 

Q ^WL[J^) q 



and the summations run over all values generated in the restricted energy 
subspace {Ei,E2). Finally, the approximate DOS used in Eq. (2) is the final 
"most accurate" DOS generated via the WL iteration [G{E) = Gwl{E)] or 
alternatively from the accumulated BH approximation [G{E) = Gbh{E)]. The 
initial modification factor of the WL process is taken to be /i = e = 2.718 .. . 
and, as usual, we follow the rule fj+i — ^JYj and a 5% flatness criterion [46,47]. 
Details on the N-fold implementation can be found in Refs. [47,49,65,66]. As 
mentioned above, the recording of the appropriate histograms takes place in 
the high levels (j = 18 — 26) of the WL process, where the incomplete detailed 
balance condition has no significant effects. 

The CrMES restriction may be defined by requesting a specified accuracy on 
a diverging specific heat (and/or on a diverging susceptibility) [47]. Alterna- 
tively, the energy density function may be used to restrict the energy space. 
This latter restriction is very simple for presentation reasons (see below) but 
we always use the original recursion method described in Ref. [46], which from 
our experience produces the most accurate and safe location of the dominant 
energy subspace. Consider a temperature of interest T£ and let E be the value 
maximizing the probability density, at this temperature. The end-points {E±) 
of the dominant energy subspaces may be located by the above mentioned 
simple energy density function condition 

B.:^<n (4) 



where r is chosen to be a small number, independent of the lattice size (i.e. 
r = 10~^). Similarly, let M be the value maximizing the order-parameter 
density at some pseudocritical temperature, for instance the pseudocritical 
temperature of the susceptibility {T^ — Tl^^^^J. Then, the end-points {M±) 
of the dominant magnetic subspaces may be located by a similar condition 

M± : -^4^ < r. (5) 



The finite-size extensions of the above dominant energy and order-parameter 
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subspaces satisfy respectively the scaling laws of the specific heat and suscep- 
tibility with exponents ajv and 7/1/, respectively [28,29,46,47]. 

The performance limitations of entropic methods, such as the WL random 
walk and the reduction of their statistical fluctuations have recently attracted 
considerable interest [52,67,68]. For the CrMES entropic scheme, presented 
above, comparative studies using various implementations and the Metropolis 
algorithm were presented in Refs. [46,47] for the d = 2 and d = 3 Ising 
models and also in Refs. [28,29] for the square SAF model. In particular, the 
effects of the used range of the WL iteration levels on the magnetic properties 
of the system and also the effect of some refinements of the WL algorithm 
were discussed. However, for large systems statistical fluctuations significantly 
increase and their reduction demands, as usually, multiple measurements, i.e. 
carrying out many WL random walks. The recordings of the {E, (5)-histograms 
in the high levels of the WL iteration process will in general produce quite 
accurate estimates for the magnetic or other properties, but it should be noted 
that this accuracy may well depend on the lattice sizes used but also on 
the particular model simulated. From our experience, all ingredients of the 
entropic schemes should be tested and treated with caution. There are, for 
instance, several ideas in the literature for reducing the modification factor of 
the WL process but also for the proper use, i.e. recording enough statistics, 
of the histogram-fiatness criterion of the WL method [47,50,52,68]. 

Taking all these into account, we first applied a multi-energy-range approach [50] 
using the one-run entropic scheme of Ref. [47], recording the {E, M) -histograms 
only in the final levels j = 18 — 26 of the WL process. For all lattice sizes 
(L = 30, 40, ... , 100, 120, 140, 160, 200, 240, 360 and 480) this simulation was 
repeated for 10 independent random walks and averages of the resTilting DOS 
and histograms were constructed. These time- demanding simulations have 
been greatly facilitated by the dominant energy subspace restriction applied. 
Even by using an r = 10~^ accuracy level, we had to simulate only a small 
part of the order of 4 — 5% of the energy spectrum for the larger lattices. For 
L = 240 we determined the DOS for about 5400 energy levels (4.7%), which is 
about 10% wider than the resulting dominant subspace of 5000 energy levels. 
For L — 480 we carried out the simulation in a subspace of 19000 energy lev- 
els (4%) and the resulting dominant subspace was about 15500 energy levels 
(3.36%). These numerical data were used to observe the first-order nature of 
the present model and in fact the resulting evidence was decisively supporting 
the first-order character of the transition (note that Figs. 1 and 2 in the next 
Section are using also these data). However, due to the first-order character 
the statistical fluctuations for the present model were rather large, compared 
to other simpler models studied earlier by the same scheme [28,46,47,48]. Ac- 
cordingly we have carried out a second-stage entropic sampling scheme only for 
the lattices L = 40, . . . , 240. In this scheme we performed two parallel random 
walks, one using the j — 26 modification factor of the WL method (WL-walk) 
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and the other using the Lee correction [39,51,52] of the DOS (Lee-walk). The 
total duration of this process was set equal to 10 x T2q, where is the time 
needed in the original scheme for the saturation of histogram fluctuations [52] 
of the WL scheme at the level j — 26. Details of these additional simulations 
are given below. 

The second-stage simulation was repeated three times for each lattice size, 
every time starting with a different DOS, selected from our first one-run ap- 
proach. In time intervals equal to the time-unit T26, the resulting DOS's from 
the two parallel walks (WL-walk and Lee-walk) were recorded, as were also 
recorded the corresponding "microcanonical estimators" for the application of 
the BH or transition matrix methods. To continue in the next time-step (of 
duration T26) we used for both parallel walks the time-accumulated average 
DOS of the Lee-walk. Of course, this practice is one of many other possi- 
bilities and ideas for "refreshing" the DOS, that could be used to continue 
the further refinement of the DOS via a WL-modification (j = 26) or via 
a Lee-adjustment [51,52]. One can even mix the DOS's of the two parallel 
walks and/or combine with the DOS's resulting from the BH or transition 
matrix methods to continue the process (such an example was considered by 
Shell et al. in Ref. [64]). At the end of the process, we have observed the 
time development of various thermal properties of the system at some conve- 
nient temperatures, for instance the specific heat and fourth-order cumulant 
pseudotransition temperatures, as well as at a temperature close to the bulk 
transition temperature (Tc = 1.80845). The inspection of their time-behavior 
verified the accuracy of our original estimates and convinced us that the above 
described practice yields a steady improvement of the time-accumulated av- 
erage DOS of all involved methods. For large lattices, the recorded during the 
WL-walk "microcanonical estimators" , for the application of the BH method, 
appeared to suffer from some "odd fiuctuations" , possibly due to a weak vio- 
lation of the detailed balance condition, that separated their various thermal 
estimates from all other estimates obtained from the other involved methods. 
Although these "odd fluctuations" were within the error-range, we have ex- 
cluded the corresponding estimates from our thermal averages and we have 
also used only the {E, M)-histograms recorded via the Lee-walk for the fur- 
ther refinement of our susceptibility estimates. The new averages were used 
as refined estimates of all properties studied and are presented in our figures 
using the range L = 40 — 240. 

Closing this Section, let us consider two different definitions for the order- 
parameter of the model under study. First, with the help of four sublattices 
of the SAF ordering one may define a two-component order-parameter and 
finally use its root- mean-square (rms) [34] 

= {Ml + M2- {Ms + M4)} /4 
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M(2) = {Ml + M4 - (M2 + M3)} /4 



(6) 



An alternative, and numerically more convenient, definition will also be consid- 
ered from the sum of the absolute values of the four sublattice magnetizations 



The resulting behavior is very similar and in particular the finite-size exten- 
sions of the resulting CrMMS completely coincide. Therefore, for large lattices 
only the second order-parameter was used in order to minimize computer mem- 
ory requirements. Using now the above definition for the order-parameter we 
may also define the magnetic susceptibility x ^ follows: 



3 General aspects of the first-order transition 

In general, first-order transitions [4] are characterized by discontinuities of an 
order-parameter, like the internal energy or the magnetization, in an idealized 
infinite system. However, in a finite-system, a true phase transition can not oc- 
cur. Instead, the jump in the order-parameter is smoothed out into a rounded 
transition which becomes increasingly sharp as the finite dimensions of the 
system go to infinity. In this case, the aim of FSS theory is to estimate the 
width and the possible shifts of the rounded transition and to determine the 
associated scaling functions. In the following subsections we summarize the 
main difficulties in making a clear distinction between first and second-order 
transitions and outline useful aspects of the existing theories of first-order 
transitions. 

3.1 On the nature of the phase transition for the triangular SAF model 

As pointed out in the introduction, Rastelli et al. [14] considered the triangular 
model with nearest- and next-nearest-neighbor antiferromagnetic interactions. 

Using a conventional (Metropolis) Monte Carlo approach, these authors con- 
cluded that the phase transition from the ordered phase at low temperatures 
to the high temperature paramagnetic phase is first-order. The identification 
of the nature of the transition was based mainly on the clear presence of the 
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M = Y,\Mi\/A. 



(7) 




(8) 
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energy double-peaks reported. Additional indications were detected from the 
behavior of the peaks of the Binder's fourth-order energy cumulant and also 
from the scaling of the specific heat, the susceptibility, and the correspond- 
ing pseudotransition temperatures. Prom the FSS behavior of the peaks of 
the Binder's fourth-order energy cumulant they estimated limiting minimum 
values different from the value of 2/3 (expected for a continuous transition) 
which are typically expected for first-order transitions. The shift behavior of 
the pseudotransition temperatures was also found in agreement with what is 
expected at a first-order transition {T* + hL~'^). However, some prob- 

lems were encountered mainly with the FSS of the specific heat peaks. The 
estimated divergences were found in marked difference with the expectations 
of the FSS theory of first-order transitions. For the case studied here (cor- 
responding to the interaction ratio R = 1) the specific heat exponent a/v 
(C* hL°'/^) was estimated to have a value ~ 1.6 instead of the value d — 2, 
which is the common expectation of all existing theories (see below). 

A preliminary qualitative comparative study was recently presented by the 
present authors [29] where the behavior of the present triangular SAF model 
was contrasted to the analogous square SAF model. In this latter study, the 
marked difference between thermally driven first- and second-order transitions 
was illustrated by using the rather traditional method of observing differ- 
ences in the energy's and order-parameter's cumulant-behavior between the 
two models. For the triangular model the behavior was indicative of a first- 
order transition, in agreement with Ref. [14]. For the square model [R = 1) 
the behavior was strongly indicating a second-order critical point [28,29], in 
disagreement with the mean- field prediction of Lopez et al. [41,42,43] of a 
first-order transition for R < 1.2. Consequently, the presence of the energy 
double-peaks and the scaling shift behavior L^'^. in the case of the triangular 
SAF model [14,29] , are indications of a first-order transition, while the absence 
of energy double-peaks and the FSS with critical exponents different from the 
lattice dimensionality, in the case of square SAF model [28,29], are indications 
of second-order phase transition. However, caution should be paid in both 
cases to ensure that the lattice size is sufficiently large so that irrelevant fields 
have scaled away and the observed behavior is the true asymptotic one and 
not a strong finite-size effect that may cease to exist in the thermodynamic 
limit. At this point we may recall some already known pitfalls. The existence 
of the energy double-peaks for the finite lattice Baxter- Wu model [48,69] is a 
counter-example since this model undergoes a second-order phase transition. 
Such continuous phase transitions with a convex dip in the microcanonical 
entropy, including the Baxter- Wu model and the q — A Potts model in two 
dimensions, have been discussed recently by Behringer and Pleimling [70]. 
In these cases the pseudosignatures of "first-order transitions" are finite-size 
effects. Similarly the apparent "first-order transition" observed in the fixed 
magnetization version of the Ising model is also a finite-size effect [71,72,73] 
and it has been shown that this "transition" ceases to exist in the thermody- 
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namic limit. 



Hence, in order to provide strong evidence and to elucidate the distinctive 
features for models undergoing first-order transitions, accurate and detailed 
numerical studies are necessary. Such studies have been carried out on the q — 
5, g = 8, and g = 10 Potts model and most of the existing theories have been 
tested and verified using this model [15,16,21,26,27]. The demonstration of the 
first-order transition features is more difficult, as should be expected, for weak 
transitions (the q — h Potts model is a well-known example) since in these 
cases the correlation length is relatively large and strong finite-size effects may 
obscure the true asymptotic behavior [16]. The present extensive numerical 
study, is an attempt to repeat some of these tests on the triangular SAF 
model, which as will be seen is an interesting weak first-order transition. Our 
accurate numerical data will permit us to adequately determine the features of 
the first-order transition of the model, to make interesting comparisons with 
previous studies, and to draw conclusions related to the applicability of the 
existing FSS theories. 



3.1.1 The Lee-Kosterlitz method and the latent heat of the transition 

We proceed in interpreting our numerical data using the free energy bar- 
rier method proposed by Lee and Kosterlitz [15,16]. With this method one 
may unambiguously identify the order of the transition and also evaluate the 
latent heat of the transition. By computing the energy distribution at vari- 
ous temperatures, we locate a pseudotransition temperature Th, at which the 
two equal-height peaks of the energy distribution are observed. These two 
peaks arc located at the values of the energies per site eo{L) = Eo{L)/L'^ and 
e(i{L) = E(i{L)/L'^, corresponding to the ordered and disordered states respec- 
tively, and are separated by a minimum at emin{L) = Emin{L)/ L'^. Following 
Lee and Kosterlitz [16], the free energy barrier is defined with the help of the 
microcanonical free energy [16,27], -F(e, L) = — In P(e, L), where P(e, L) is the 
energy distribution at the pseudotransition temperature Th- At a first-order 
transition the microcanonical free energy is assumed to have an expansion of 
the form 



where /(e) is the bulk free energy which is a minimum and constant for to < 
e < Sd and /^(e) is a surface free energy which has a maximum at emin- The 
arguments presented by Lee and Kosterlitz [16] suggest the following scaling 
forms for the finite-size estimates of the locations of the equal-height peaks 



F(e,L) = L'^/(e) + L'^-V.(e) + ---, 



(9) 



(10) 
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and 



ediL) = ed + 0{L-'), 



(11) 



where Cq and are the bulk energies of the ordered and disordered states 
respectively. From the above assumptions for the bulk and surface free ener- 
gies, it follows that the maximum of the surface energy at Cmin produces a free 
energy barrier that scales according to 



This scaling form has been used as a basic test [15,16,27,74] for detecting a 
temperature-driven first-order transition, while the predicted scaling behavior 
[Eqs. (10) and (11)] can be used to extrapolate and obtain the latent heat of 
the transition Ae — 64 — 60. 

Fig. 1 presents the energy distributions for small and large lattice sizes. Al- 
though, for the larger lattices the simulation error fluctuations are evident, it 
is also very clear that with increasing lattice size the barrier between the two 
peaks is steadily increasing, signaling the emergence of the expected two delta- 
peak behavior in the thermodynamic limit. Thus, this figure is an effective 
demonstration of the persistence of the first-order transition. Fig. 2 (a) is the 
corresponding illustration, including all lattice sizes studied, of the behavior 
of the free energy barrier, defined above. It is apparent from this figure that 
the surface tension, that is the quantity AF/L — {kBT\a{Pjnax/ Pmin)}/L, 
tends to a non-zero value and the attempt to extrapolate by the dotted line 
gives a value of the order 0.00539(12). Thus, the expectation (12) of the ar- 
guments of by Lee and Kosterlitz [16] is well verified. Fig. 2(b) illustrates 
the L-dependence of the location of the energy density peaks (minima of the 
free energy) as well as their difference. Again in the scale shown, the expec- 
tations described by Eqs. (10) and (11) are well satisfied. Thus, using linear 
extrapolations in the range of large L, we have estimated the bulk energies 
60 = -1.61784(728) and = -1.48852(565), and the latent heat of the transi- 
tion Ae = — Co = 0.129(12). This is to be compared with the value 0.176(6) 
estimated by Rastelli et al. [14]. 

3.2 The double Gaussian approximation 

K. Binder and co-workers [20,21,75] introduced the so-called double Gaus- 
sian approximation. This assumption concerns the probability distribution of 
the order-parameter and has become a widely used feature in the theories of 
first-order transitions. The original approach leading to this approximation 



AF(L) = F{e^in, L) - F{eo, L) - L 



d-l 



(12) 
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was based on the intuitive extension of the Gaussian approximation of a sin- 
gle phase, which may be theoretically justified by appealing to the central 
limit theorem. It was argued [20,21,75] that, the distinctive feature of a first- 
order transition is the coexistence of phases at the transition point and in fact 
this important observation is the starting point of other theories of first-order 
transitions [15,16,24,20,21,26,27,75,76,77,78,79,80,81,82,83,84,85]. Ignoring ef- 
fectively interface effects between domains of different phases, the equilibrium 
of the phases is described by a superposition of Gaussians centered at different 
values of the order parameter corresponding to the different coexisting phases. 

Furthermore, it was soon realized [15,16] and theoretically [78,79,80,81,82,83,84,85] 
justified, that the same qualitative picture can be obtained by a more fun- 
damental assumption. The partition function of a large finite system, with 
periodic boundary conditions, undergoing in the limit L — > oo a first-order 
transition, is well approximated at the transition point as a sum of terms, 
each describing a separate phase. For temperature-driven first-order phase 
transitions, in which q equivalent ordered phases coexist at the transition 
temperature (/3c) with one disordered phase, the dominant contributions in 
the partition function, separated in contributions from the different coexist- 
ing phases, may be written as [16] 

Z ^ ^g-/3c/i(/3c)i'*^ 
1=1 

where /i(/3c) is the free energy per site of the ith phase in the thermodynamic 
limit. This fundamental relationship is a statement that each phase contributes 
equally to Z at the transition temperature {fi{Pc) — fiPc))- As pointed out by 
Lee and Kosterlitz [16], the double Gaussian approximation may be obtained 
by assuming an extension of (13) in the vicinity of Tc and using appropriate 
temperature Taylor expansions for the free energies /j(/3) of the coexisting 
phases. The probability density for the energy can then be obtained by the 
inverse Laplace transform of the expanded partition function and, in order 
0{L'~'^) in the free energy expansion, one can derive the double Gaussian 
form 

where $G(e; Cji (^i) is the normal distribution with mean Cj = ej+Cj(^)(T— Tc) 
and variance = ksT'^CiL''^ . and Cj are the energy and specific heat of 
the ith bulk phase, derived as appropriate derivatives from the corresponding 
free energies. Furthermore, 

A = 1 Ing + i(e„ - e,)^{T - T^)L'' + \{C^ - C,) (^y^)' (15) 



(13) 
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so that the ratio of weights corresponding to the two terms of Eq. (14) satisfies 
the relation 

= e^^ = ^e^(/.-/o)i^ (16) 



From the above distribution, or in a more straightforward way, by using dif- 
ferentiation of the expanded partition function one can deduce the finite-size 
behavior of all energy cumulants, as shown by Lee and Kosterlitz [16]. These 
tedious calculations have been also repeated by Janke [27] and in the follow- 
ing we summarize only the relevant results that have been used in our study 
of the triangular SAF model. At this point let us emphasize the fact, stated 
explicitly bellow, that the double Gaussian approximations predicts, in gen- 
eral, that all finite-size contributions enter in the scaling equations in powers 
of L"*^. The general shift behavior predicted for the various pseudotransition 
temperatures (denoted collectively as T^*) is 

T* = T^ + bL-'^{l + 0{L-'')), (17) 



where these temperatures correspond to the peaks of various energy cumu- 
lants such as those defined bellow. As noted above, not only the main shift 
contribution but also the higher-order corrections enter in Eq. (17) only in 
powers of L'^. Three different energy cumulants have been used in our calcu- 
lations, the specific heat per site C(L), the second-order cumulant U2{L), and 
the Binder's fourth-order cumulant V4(L). Bellow we give their definitions and 
the predicted scaling behavior [16,27] at the corresponding pseudotransition 
temperature, at which the cumulant peak (maximum or minimum) occurs 

C{L) = L-'' [{E^) - {Ef] /T^ (18) 
C{L)L-'^U^ = + 0{L-'j) (19) 

c 

U^iL) = (20) 

U2{L)\ma. ^ ^-^^-^ + 0{L-') (21) 

V4L)Un + 0{L-% (23) 
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In the following Section we proceed to test the above scaling laws for the trian- 
gular SAF model. In particular, the coefficient of the dominant diverging term 
of the specific heat and the limiting values of the cumulants U2 and V4, will be 
determined in two independent ways: (a) by finding their values from convinc- 
ing fitting assumptions of our peak-data and (b) from the extrapolated values 
of the energies Co and (see Fig. 2(b) and Eqs. (10) and (11)) and the above 
predictions. The proposed comparison will provide a useful self-consistency 
test of our results but also of the theoretical predictions. It will also assist 
a critical discussion on the higher-order corrections of the double Gaussian 
approximation. Noteworthy, that the above leading terms may also be derived 
heuristically from a simple two- phase model, as shown by Janke [27] . 

3.3 Number- of -phases parameters and the determination of transition points 

Traditionally, the general shift behavior T* k, Tc-\-bL~'^ has been used for both 
the identification of a first-order transition and the determination of the tran- 
sition temperature. However, Borgs and Janke [26] have suggested additional 
methods facilitating the determination of transition points for systems with 
periodic boundary conditions. Two such methods will be discussed below and 
then elaborated in the sequel to observe the behavior of the present first-order 
transition. In the first method, a parameter is introduced that requires two 
different finite lattices of volumes Vi — Lf and V2 — respectively. This 
parameter has been called the number-of-phases parameter [27] and is defined 
as: 

where a is the volume ratio a — V2/V1 — {L2/L1Y. The origin of this param- 
eter is the fundamental relationship (13), which is extended, as earlier, in the 
vicinity of the transition point by introducing [26,27] metastable free energy 
densities /i(/5). These are defined in such a way that are equal to the bulk free 
energy /(/3) in the temperature range in which the corresponding phase is a 
stable phase and strictly larger than /(/3) when the phase is unstable. The 
coexistence partition function is now written as [26,27] 

Z,er{V; P) = |e e-^^^(^)^ I (1 + 0(ye-^/^i) . (25) 

A heuristic derivation of the above exponential correction bound has been 
presented in Ref. [26]. The hypothesis for exponential small shifts, to be dis- 
cussed in the next Section, derives from the above bound. Eq. (25) and the 



Zper{V2;P)\ 



(24) 
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assumptions of the behavior of metastable free energy densities may be used 
to show [27] that the parameter 

7V(/?)= hm Z,e.(l^;/3)e^^(^)^ (26) 



is equal to the number of stable phases, which with increasing temperature 
takes the values g, q + 1, and 1. A finite-size approximation of the above 
parameter will be expected to have a peak at a finite-size transition point 
(Tv/y) and the number-of- phases parameter defined in Eq. (24) is a suitable 
such approximation that may be used to locate the transition point. The 
condition for the maximum of the parameter of Eq. (24) corresponds to the 
crossing point of the mean energy functions of the pair of lattices involved and 
reads as: 

a{Ei)\Ty,y ^ {E2)\Ty/y or (ei)|T^ = (e2)|rw^. (27) 



The location of the peaks of the parameter (24) can be estimated by a dom- 
inant CrMES scheme, such as that followed by the present authors in other 
applications [28,46,47]. Of course, the restricted energy space scheme should 
be sufficient for both lattices of the pair in the temperature range of interest, 
that is in the region of the peaks. We shall use the sequence of pairs of lattices: 
(40, 80), (50, 100), (60, 120), (70, 140), (80, 160), (100, 200), and (120, 240), 
with volume ratio a — A (L2 = 2Li, d — 2). Fig. 3(a) illustrates the behavior 
which is similar with the behavior presented by Janke [27] for the Potts model. 

We now discuss a variant of the above number-of-phases parameter. First, we 
define a "reduced partition function" . The maximum term, at each tempera- 
ture, has been factorized and dropped out from the partition function's sum. 
Thus, 

Z{V;P)= Y: (28) 

(E-,E+) 



where E — E{T) is the most probable energy corresponding to the maximum 
term, and 

^E) = S{E)-f5E-[S{E)-(5E]. (29) 



Since the most probable energy has two values at the equal-height temperature 
(T/i), the above reduction introduces in the above definition microcanonical 
features of the transition. These features are now carried over to the new ratio 
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parameter 



A 



Nx{V,,V2;P) 



(30) 



Let us regard the exponent A as a free parameter, not necessarily equal to 
the volume ratio of the two systems (a = 4), and observe the resulting be- 
havior. Fig. 3(b) shows the behavior of this new ratio parameter for the case 
A = 2 and the three larger lattice size pairs of the sequence considered above. 
It can be seen from this figure that by increasing the temperature this pa- 
rameter shows a first maximum peak (Tj^axi) which resembles the peak of 
the original number-of-phases parameter. Further increase in the temperature 
yields, for each pair, two additional very sharp peaks, one minimum (Tmin) 
and one maximum {Tmax2)- These two peaks are reflections of the maximum 
term reduction in the definition (30). The sharp minimum corresponds to the 
equal-height pseudotransition temperature of the larger lattice, while the 
sharp maximum corresponds to the equal-height pseudotransition tempera- 
ture Tfi of the smaller lattice. The value of the exponent A determines the 
location of the first maximum peak and the sharpness of these graphs. Thus, 
a similar but less pronounced picture is obtained for the values A = 3 and 
A = 4, and it is found that the location of the two new sharp peaks are not 
notably influenced by the change in the value of the exponent A. This can be 
easily understood since the new sharp peaks are in effect properties of the two 
separated lattices. The peaks of this new ratio can be also be described by an 
approximate crossing condition, similar to Eq. (27), which may be written as: 



In the above condition we have neglected the temperature variation of the most 
probable energy value E. This variation will produce vanishing contributions, 
as tested also numerically, involving the difference between microcanonical 
and canonical temperature [{'&Si/'i!}Ei)\^ — /?]. 

Fig. 4(a) is an illustration of the crossing point of the mean energy functions of 
a pair of systems (Li = 80, L2 = 160) corresponding to the peak of the original 
number-of-phases parameter (24). In this the behavior of the most probable 
energies of the pair of systems is also illustrated. Fig. 4(b) presents the tem- 
perature behavior of the differences, (e) — e, of the two characteristic energies 
for each system of the chosen pair. For the larger lattice the graph correspond- 
ing to the double of the difference (e) — e is also shown in order to facilitate 
the illustration of the above crossing relationship (31) for A = 2. Notewor- 
thy, that this illustration provides an alternative practical way of locating the 
equal-height pseudotransition temperature Th, since at this temperature (in- 
version point) the corresponding difference changes sign. As we increase the 




(31) 
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temperature, the difference cTirve of the smaller lattice is crossed by the curve 
of the larger lattice at a point {T„iaxi) corresponding to the first maximum in 
Fig. 3(b). As can be seen the location of this point (Tx) depends on the value 
of the exponent A and for the value X — a/2 — 2, the crossing is deeper and 
occurs at a lower temperature, than in the case corresponding to A = a = 4. 
This explains why the first maximum peak (Tmaxi) in the case A = 2 is more 
pronounced as mentioned earlier. Finally, the sharp minimum {Tmin) and the 
sharp maximum {Tjnax2) in Fig- 3(b) are described by the crossing points cor- 
responding to the inversion points of the larger and smaller lattice respectively 
and their location is not practically influenced by the value of the exponent A. 

We shall conclude this Section with a brief reference on the second definition 
of a finite-size transition point introduced by Borgs and Janke [26] . This equal- 
weight finite-size transition point (Tw) requires data from one lattice only and 
is the point where the ratio of the total weight of the ordered phases to the 
weight of the disordered phase approaches q, so that 

R{V,E)= P{E)/ P{E) ^"^W^ = (32) 

E<Em E>Em ^^<^ 

where Ej^ may be chosen in various ways and following Refs. [26,27] we will 
choose it to be the energy of the minimum between the two equal-height peaks 
at the pseudotransition temperature Th- As a variant of the above definition, 
we have attempted to satisfy Eq. (32) using the part of the weights above the 
minimum between two unequal peaks. The resulting transition points have 
similar behavior and will not be presented. 



4 Analysis of numerical data 

4.1 Estimation of the transition temperature 

Using our accurate data and various alternative methods we have estimated 
the transition temperature with an error at most in the fifth significant fig- 
ure. Our safe estimate is = 1.8085(1) and is an improvement of the three 
estimates of the order of Tc = 1.8078(1), given originally by Rastelli et al. [14] 
and close to their magnetic fourth-order cumulant estimate = 1.8084(1). 
An interface method estimation of this transition point {Tc — 2.044) was given 
by Slotte and Hemmer [86]. Fig. 5(a) presents a first attempt to estimate the 
transition temperature by a simultaneous fitting on five pseudotransition tem- 
peratures. As indicated, in an obvious notation on the graph, these are the 
equal-height temperature (T/j) and the pseudotransition temperatures corre- 
sponding to the peaks of the specific heat, second- and fourth-order energy 
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cumulants, and finally the peak of the snsceptibility. In the first attempt the 
shift behavior predicted by the double Gaussian approximation has been as- 
sumed including also the first shift correction term (i.e. 0{L~'^)). Fig. 5(b) 
presents a second attempt using now a variant of the above assumption by 
replacing the shift correction term with an correction. This replacement 
makes practically no difference in estimating the transition temperature be- 
cause the dominant shift term is about three orders of magnitudes larger than 
the next correction in all cases, and the asymptotic behavior has been well 
reached in the temperature shift behavior. Noteworthy that, if we restrict 
the fitting attempt in the range L = 30 — 140, we will obtain an estimate 
Tc = 1.80845(5). 

Next we present the behavior of the new transition point observables intro- 
duced by Borgs and Janke [26], as well as the similar ones introduced and 
discussed in the last subsection. Although we will not try to further refine 
the above accurate estimate, it is of interest to illustrate their behavior, show 
that they lead to the same estimate, and discuss their behavior in comparison 
with previous studies on the Potts model. Fig. 6(a) illustrates the behavior of 
the three peaks of the reduced number-of-phases parameter defined in the last 
Section in Eq. (30) for A = 2. The data for the two peaks {T^ax2 and T^i„) are 
approaching from above the transition point and coincide in practice with the 
equal-height temperatures for the smaller and larger lattices respectively, as 
explained in subsection 3.3. The data for the peaks Tmaxi, which resemble the 
original peaks of the number-of-phases parameter of Borgs and Janke [26] , are 
also moving steadily to the expected limit and their fitting, using Eq. (17), 
yields again the estimate = 1.8085(1). Fig. 6(b) collects the behavior of 
six finite-size transition points and compares them with our estimate for the 
transition point and also with the estimate Tc = 1.8078 of Rastelli et al. [14]. 
We have included in this figure two transition points depending only on one 
lattice, namely the equal- height transition point (Th) and the weight transi- 
tion point (TV) defined in Eq. (32). The rest of the shown transition points 
are the original number-of-phases parameter transition point introduced by 
Borgs and Janke [26], denoted as Ty/v, and the left maximum point {Tmaxi 
or Tx) of the reduced number-of-phases parameter introduced in this paper 
for three values of the exponent A, A = 2,3 and 4 (see the discussion below 
Eq. (31)). Their limiting behavior is an extra verification of our estimate for 
the transition temperature. 

Borgs and Janke [26] in their Fig. 2, and Janke [27] in his Fig. 11, have 
illustrated the behavior of the transition points Ty/v and the 3V for the Potts 
model with q — 5, q — 8, and q — 10. As it can be seen from the mentioned 
figures, the q = 8 and q = 10 finite-size transition points approach their bulk 
limit for relatively small lattices. Therefore, these authors have concluded 
that for strong first-order phase transitions the exponential bound in Eq. (25) 
provides a natural explanation of the observed "exponentially small" shifts 
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with respect to the infinite- vohime transition point. For the weak first-order 
transition of the g = 5 Potts model, the weight transition point T^^/ can also 
be seen from the mentioned figures to have reached the bulk limit even for 
the smallest lattice size L — 2d. However, the number-of-phases parameter 
transition points Ty/v did not obey this exceptional behavior, but rather its 
shift behavior resembled the shift behavior of the traditional pseudotransition 
temperatures, i.e. an L~'^ shift. Accordingly, Janke [27] has raised also the 
question of a possible fortuitous cancellation and has suggested the interest 
to the investigation of these exponentially small shifts in other weak first- 
order transitions. Fig. 6(b) provides an example showing that for the present 
weak first-order transition the expectation for exponentially small shifts is 
not verified for both transition points Tw and Ty/v- Their behavior differs 
appreciably from their bulk limit and the reduced number-of-phases parameter 
for the value A = 4 yields the smaller shifts, smaller than that of the natural 
number-of-phases parameter. 



4-2 Behavior of the energy cumulants and the magnetic susceptibility 

It is well known that the FSS behavior of the specific heat, and in particu- 
lar its peak-scaling behavior, may be obscured by the presence of unknown 
correction terms and the extraction of the scaling exponent may be a numer- 
ically indecisive and difficult task [39,40,74]. The accuracy of numerical data 
is then of crucial importance [40] . This unpleasant problem may also occur in 
the estimation of other related energy cumulants. For the present model such 
problems were obviously encountered by Rastelli et al. [14]. Their estimation 
of the specific heat exponent a/u {C* ~ bL°'^'^) was found in clear disagree- 
ment (estimated to have a value ~ 1.6) with the theoretically expected value 
a/v — d — 2. Since the evidence for the present first-order transition is now 
overwhelming we will adopt the generally accepted theoretical prediction for 
the leading contribution to the specific heat divergence. Thus, we shall fix the 
value of the exponent a/ v = d = 2 and we will examine whether there is a 
stable form of corrections terms that is compatible with our accurate data in 
the peak-region. 

Fitting our numerical data from the first 10 WL runs in various L-ranges, 
we found that the simple correction term of the form cL~^, instead of cL~^ 
predicted by the double Gaussian approximation, yielded stable fittings in all 
the L-ranges tried, from L = 30 to L = 480. The very small variations in 
the involved coefficients (5% or less) of the fitting attempts are, of course, 
strong indications in favor of this otherwise at hoc assumption. Fig. 7 shows 
that even the remote L-ranges: L = 30 — 140 and L = 100 — 480 produce 
4% deviation in the coefficient of the leading term and only 1% in the coeffi- 
cient of the assumed correction term. Fig. 8 provides further strong support 
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in favor of the above assumption using our second-stage refined data in the 
range L = 40 — 240. This figure present simultaneous fittings of the specific 
heat data not only at its pseudotransition temperature but also in the other 
pseudotransition temperatures corresponding to the peaks of the other energy 
cumulants, the peak of the magnetic susceptibility and the equal-height transi- 
tion point. Noteworthy, that such a simultaneous fitting approach, quite easily 
implemented within our entropic scheme, would be a rather demanding task 
in a traditional importance sampling Monte Carlo scheme. In this way we are 
probing systematically the FSS behavior of the divergence of the specific heat 
in a wide region around its peak. Fig. 8(a) illustrates the fitting attempt based 
on the above mentioned (cL~^) assumption and Fig. 8(b) shows a remarkable 
decline from the double Gaussian prediction (cL~^). The double Gaussian pre- 
diction does not correctly describe the correction term and this is reflected in 
the resulting very large error value which is about 35 times larger than the 
value corresponding to the well obeyed correction term {cL'^). Let us now 
use the corresponding coefficients of the dominant divergences for the peaks 
of the specific heat to compute the resulting values for the latent heat and 
compare these values with the value obtained earlier from the extrapolations 
of Fig. 2(b), based on the proposal of Lee and Kosterlitz [16] and described in 
Eqs. (10) and (11). The resulting values for the latent heat from the fittings 
in Fig. 8 are in the range 0.121 — 0.133 and 0.143 — 0.163 respectively. These 
values should be compared with the value 0.129(12), see also Table 1 below, 
obtained from the extrapolated values of Cq and e^. 

An analogous behavior is observed in the FSS cbehavior of the other two 
energy cumulants. Fig. 9(a) shows in the same graph fitting attempts using 
the data for the minima of the Binder's fourth-order cumulant presuming 
again the two different assumptions for the correction terms (6L^^ and hL~^). 
As can be seen from this figure the double Gaussian prediction fails again to 
correctly describe the correction term. The resulting very large error value 
is now 135 times larger than the value corresponding to the simple correction 
term of the form hL^^ . Furthermore, the limiting value corresponding to the 
assumption using the double Gaussian correction term is 0.66259(40) for the 
L-range shown, which differs appreciably from the value 0.66435(33), obtained 
from Eq. (23) and the extrapolated values of to and td given in Table 1. On the 
other hand the limiting value resulting from a leading correction of the form 
gives 0.66465(25) for the L-range shown, as illustrated in Fig. 9(b), and 
favors an expansion starting with the term hL^^ . Fig. 10 repeats the above 
comparison for the second-order cumulant U2- The expansion starting with 
the term hL"^ is surprisingly well obeyed, while the double Gaussian behavior 
is inconsistent. The limiting value obtained from double Gaussian correction 
term is 1.00308(25) in clear disagreement with the value 1.00174(25) obtained 
from Eq. (21) and the extrapolated values of Cq and e^, while the proposed 
expansion in Fig. 10(b) gives the value 1.00163(15) quite in agreement with 
the also approximate estimate 1.00174(25). 
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Table 1 summarizes the estimates for the above thermal characteristics includ- 
ing also the corresponding exact values [16,27,87,88] for the cases q = 10, 8, 6 
and g = 5 of the Potts model. The estimates for the triangular SAP model 
with R—1 interpolate between the values of the q — 5 and q — 6 Potts model. 
Therefore, we have called the present transition a weak first-order transition. 
It appears that the six- fold degeneracy of the ground state of the present 
model plays an important role in determining the weakness of the transition, 
and possibly the occurrence of its first-order character. Pinally, let us turn 
to the behavior of the magnetic susceptibihty x- We present in Pig. 11 the 
divergence of the susceptibihty x at the five pseudocritical temperatures T*, 
using again the two types of corrections, i.e. an (Pig- 11(a)) and an 
(Pig. 11(b)), respectively. The solid lines are simultaneous fitting attempts, 
using the power laws shown in the figures. Comparing Pigs. 11(a) and (b) we 
see that, although one can not just by inspection differentiate between an 
and an correction term, is smaller for the former correction. 



5 Conclusions 

The triangular Ising model with nearest- and next-nearest-neighbor antiferro- 
magnetic interactions has been studied and its first-order transition features 
have been clarified, when the interaction ratio is i? = 1. We have outlined the 
most important aspects of the existing theories of first-order transitions and 
we have tested on this model some basic hypothesis from these theories, com- 
paring our results and findings with previous work on the Potts model. Our 
numerical data have been used to obtain accurate estimates for all the ther- 
mal characteristics of the present weak first-order transition. All the generally 
accepted predictions of the finite-size scaling theory for first-order transitions, 
concerning transition point shifts, thermal, and magnetic anomalies, have been 
well verified for the present model. 

However, two of our findings for this model are not compatible with some 
theoretical or phenomenological expectations. The first of these concerns the 
behavior of the new transition point observables introduced by Borgs and 
Janke [26] and also some similar transition points introduced in this paper. 
These finite-size transition points are suitable finite-size approximations of 
the fundamental number-of-phases parameter and it has been theoretically 
argued [16,27] that should be expected to obey exponentially small shifts. 
This expectation is not verified for the present weak first-order transition. 

Pinally, we have shown that the well-known double Gaussian approximation 
does not describe correctly the higher correction terms for all energy cumu- 
lants of the present model. It appears that the first correction term in the 
expansions of energy cumulants is of the form L~'^+^(= L~^) and not of the 
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form L~'^{= L~'^), expected from the double Gaussian approximation. Lee and 
Kosterlitz [16] have pointed out the inadequacy of the Gaussian approxima- 
tion to produce shifts in the locations of the energies of the equal-height peaks 
in agreement with those described by Eqs. (10) and (11), which are rather well 
observed in simulations (see Ref. [16] and Fig. 2(b)). In fact the shortcomings 
of the Gaussian behavior has been critically discussed in the early work of 
Challa et al. [21]. Our analysis is therefore suggesting again the need for a 
more realistic theory. It is tempting to assume that the attempted here and 
well obeyed expansions for the energy cumulants, starting with the correction 
term bL"^, may have their origin in the neglected domain wall corrections 
and could have the same explanation with the existing puzzling situation con- 
cerning the shift behavior of the free energy minima pointed out by Lee and 
Kosterlitz [16]. 



Acknowledgements 

This research was STipported by the Special Account for Research Grants of 
the University of Athens under Grant Nos. 70/4/4071. N.G. Fytas would like 
to thank the the Alexander S. Onassis Public Benefit Foundation for financial 
support. 



References 

[1] W. Selke, in Phase Transitions and Critical Phenomena, edited by C. Domb and 
J.L. Lebowitz, Academic Press, London, 1992, vol. 15. 

[2] D. Lawrie, S. Sarbach, in Phase Transitions and Critical Phenomena, edited by 
C. Domb and J.L. Lebowitz, Academic Press, London, 1992, vol. 9. 

[3] L. Onsagcr, Phys. Rev. 65 (1944) 117; R.J. Baxter, Exactly solved models in 
statistical mechanics, Academic Press, London, 1982; B. McCoy, T. Wu, The 
Two Dimensional Ising Model, Harvard University Press, Cambridge, 1972. 

[4] K. Binder, Rep. Prog. Phys. 50 (1987) 783. 

[5] Y. Tanaka, N. Uryii, J. Phys. Soc. Jpn. 39 (1975) 825. 

[6] M.K. Phani, J.L. Lebowitz, M.H. Kalos, Phys. Rev. B 21 (1980) 4027. 

[7] J. Oitmaa, J. Phys. A 14 (1981) 1159. 

[8] M.J. Velgakis, J. Oitmaa, J. Phys. A 21 (1988) 547. 

[9] D.P. Landau, Phys. Rev. B 27 (1983) 5604. 



24 



[10] M. Bretz, Phys. Rev. Lett. 38 (1977) 501. 

[11] E. Domany, M. Schick, J.S. Walker, R.B. Griffiths, Phys. Rev. B 18 (1978) 2209. 

[12] H.H. Roomany, H.W. Wyld, Phys. Rev. B 23 (1981) 1357. 

[13] P. Peczak, D.P. Landau, Phys. Rev. B 39 (1989) 11932. 

[14] E. Rastelh, S. Regina, A. Tassi, Phys. Rev. B 71 (2005) 174406. 

[15] J. Lee, J.M. Kosterhtz, Phys. Rev. Lett. 65 (1990) 137. 

[16] J. Lee, J.M. Kosterhtz, Phys. Rev. B 43 (1991) 3265. 

[17] For a review of FSS see M. N. Barber, in Phase Transitions and Critical 
Phenomena, edited by C. Domb and J.L. Lebowitz, Academic, New York, 1983, 
Vol. 8. See also K. Binder, ibid., Vol. 8. 

[18] M.E. Fisher, in Critical phenomena, Proc. Enrico Fermi International School of 
Physics, ed. M.S. Green, Academic Press, New York, 1971, Vol. 51. 

[19] M.E. Fisher, M.N. Barber, Phys. Rev. Lett. 28 (1972) 1516. 

[20] K. Binder, D.P. Landau, Phys. Rev. B 30 (1984) 1477. 

[21] M.S.S. Challa, D.P. Landau, K. Binder, Phys. Rev. B 34 (1986) 1841. 

[22] V. Privman, in Finite-size scaling and numerical simulation of statistical 
systems, ed. V. Privman, World Scientific, Singapore, 1990. 

[23] K. Binder, in Computational methods in field theory, eds. H. Gausterer, C.B. 
Lang, Springer, Berlin, 1992, p. 59. 

[24] M.E. Fisher, A.N. Berker, Phys. Rev. B 26 (1982) 2507. 

[25] B. Nienhuis, A.N. Berker, E.K. Riedel, M. Schick, Phys. Rev. Lett. 43 (1979) 
737. 

[26] C. Borgs, W. Janke, Phys. Rev. Lett. 68 (1992) 1738. 

[27] W. Janke, Phys. Rev. B 47 (1993) 14757. 

[28] A. Malakis, P. Kalozoumis, N. Tyraskis, Eur. Phys. J. B 50 (2006) 63. 

[29] A. Malakis, P. Kalozoumis, N.G. Fytas, Rev. Adv. Mater. Sci. (in press). 

[30] D.P. Landau, K. Binder, Phys. Rev. B 31 (1985) 5946. 

[31] K. Minami, M. Suzuki, J. Phys. A 27 (1994) 7301. 

[32] M.P. Nightingale, Phys. Lett. A 59 (1977) 486. 

[33] R.H. Swendsen, S. Krinsky, Phys. Rev. Lett. 43 (1979) 177. 

[34] K. Binder, D.P. Landau, Phys. Rev. B 21 (1980) 1941. 

[35] D.P. Landau, Phys. Rev. B 21 (1980) 1285. 



25 



[36] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, Teller, J. 
Chem. Phys. 21 (1953) 1087. 

[37] A.B. Bortz, M.H. Kales, J.L. Lebowitz, J. Comput. Phys. 17 (1975) 10. 

[38] K. Binder, Rep. Prog. Phys. 60 (1997) 487. 

[39] M.E.J. Newman, G.T. Barkema, Monte Carlo Methods in Statistical Physics, 
Clarendon Press, Oxford, 1999. 

[40] D.P. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical 
Physics, Cambridge University Press, Cambridge, 2000. 

[41] J.L. Moran-Lopez, F. Aguilera-Granja, J.M. Sanchez, Phys. Rev. B 48 (1993) 
3519. 

[42] J.L. Moran-Lopez, F. Aguilera-Granja, J.M. Sanchez, J. Phys. C 6 (1994) 9759. 

[43] E. Lopez-Sandoval, J.L. Moran-Lopez, F. Aguilera-Granja, Solid State 
Commun. 112 (1999) 437. 

[44] C. Buzano, M. Pretti, Phys. Rev. B 56 (1997) 636. 

[45] R.J. Baxter, Ann. Phys. (N.Y.) 70 (1972) 193. 

[46] A. Malakis, A. Peratzakis, N.G. Fytas, Phys. Rev. E 70 (2004) 066128. 

[47] A. Malakis, S.S. Martinos, LA. Hadjiagapiou, N.G. Fytas, P. Kalozoumis, Phys. 
Rev. E 72 (2005) 066120. 

[48] S.S. Martinos, A. Malakis, LA. Hadjiagapiou, Physica A 352 (2005) 447. 

[49] A. Malakis, N.G. Fytas, Phys. Rev. E 73 (2006) 016109. 

[50] F. Wang, D.P. Landau, Phys. Rev. Lett. 86 (2001) 2050; Phys. Rev. E 64 (2001) 
056101; B.J. Schulz, K. Binder, M. Muller, D.P. Landau, ibid. 67 (2003) 067102. 

[51] J. Lee, Phys. Rev. Lett. 71 (1993) 211. 

[52] H.K. Lee, Y. Okabe, D.P. Landau, Comput. Phys. Commun. 175 (2006) 36. 

[53] P.M.C. de Oliveira, T.J.P. Penna, H.J. Herrmann, Braz. J. Phys. 26 (1996) 677. 

[54] J.-S. Wang, T.K. Tay, R.H. Swendsen, Phys. Rev. Lett. 82 (1999) 476; J.-S. 
Wang, R.H. Swendsen, J. Stat. Phys. 106 (2002) 245. 

[55] B.A. Berg, T. Neuhaus, Phys. Lett. B 276 (1991) 249; Phys. Rev. Lett. 68 (1992) 
9. 

[56] G.R. Smith, A.D. Bruce, J. Phys. A 28 (1995) 6623. 

[57] G.M. Torrie, J.-P. Valleau, J. Comput. Phys. 23 (1997) 187. 

[58] R.H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 57 (1986) 2607. 



26 



[59] C.J. Geyer, in Computing Science and Statistics: Proceedings of the 23rd 
Symposium on the interface, ed. E.K. Keramidas, Interface Foundation, Fairfax 
Station, New York, 1991, p. 156. 

[60] E. Marinari, G. Parisi, Europhys. Lett. 19 (1992) 451. 

[61] A. P. Lyubartsev, A. A. Martsinovskii, S.V. Shevkunov, P.N. Vorontsov- 
Velyaminov, J. Chem. Phys. 96 (1992) 1776. 

[62] K. Hukushima, K. Nemoto, J. Phys. Soc. Jpn. 65 (1996) 1604. 

[63] E. Marinari, G. Parisi, J. Ruiz-Lorenzo, in Spin Glasses and Random Fields, ed. 
A. P. Young, Directions in Condensed Matter Physics, World Scientific, Singapore, 
1998, Vol. 12. 

[64] C. Yamaguchi, Y. Okabe, J. Phys. A 34 (2001) 8781; P.N. Vorontsov- 
Velyaminov, N.A. Volkov, A.A. Yurchenko, J. Phys. A 37 (2004) 1573; N. Rathore, 
J.J. de Pablo, J. Chem. Phys. 116 (2002) 7225; M.S. Shell, P.G. Debenedetti, A.Z. 
Panagiotopoulos, Phys. Rev. E 66 (2003) 056703; P. Poulain, F. Calvo, R. Antoinc, 
M. Broyer, Ph. Dugourd, ibid. 73 (2006) 056704; C. Zhou, T.C. Schulthess, S. 
Torbriigge, D.P. Landau, Phys. Rev. Lett. 96 (2006) 120201. 

[65] B.J. Schulz, K. Binder, M. Muller, Int. J. Mod. Phys. C 13 (2002) 477. 

[66] A. Malakis, S.S. Martinos, LA. Hadjiagapiou, A.S. Peratzakis, Int. J. Mod. 
Phys. C 15 (2004) 729. 

[67] P. Dayal, S. Trebst, S. Wessel, D. Wiirtz, M. Troyer, S. Sabhapandit, S.N. 
Coppersmith, Phys. Rev. Lett. 92 (2004) 097201. 

[68] C. Zhou, R.N. Bhatt, Phys. Rev. E 72 (2005) 025701(R). 

[69] N. Schreider, J. Adler, J. Phys. A 38 (2005) 7253. 

[70] H. Behringer, M. Pleimhng, Phys. Rev. E 74 (2006) 011108. 

[71] M. Plcimling, A. Hiiller, J. Stat. Phys. 104 (2001) 971. 

[72] K. Binder, Physica A 319 (2003) 99. 

[73] S.S. Martinos, A. Malakis, LA. Hadjiagapiou, Physica A 366 (2006) 273. 
[74] D. Ledue, D.P. Landau, J. Teillet, J. Appl. Phys. 83 (1998) 6305. 
[75] K. Binder, Z. Phys. B 43 (1981) 119. 

[76] H.W.J. Blote, M.P. Nightingale, Physica A 112 (1982) 405. 
[77] V. Privman, M.E. Fisher, J. Stat. Phys. 33 (1983) 385. 

[78] C. Borgs, R. Kotecky, J. Stat. Phys. 61 (1990) 79; Phys. Rev. Lett. 68 (1992) 
1734; Physica A 194 (1993) 128. 

[79] C. Borgs, R. Kotecky, S. Miracle-Sole, J. Stat. Phys. 62 (1991) 529. 

[80] C. Borgs, S. Kappler, Phys. Lett. A 171 (1992) 37. 



27 



[81] C. Borgs, J.Z. Imbrie, Commun. Math. Phys. 145 (1992) 235. 
[82] C. Borgs, Nucl. Phys. B 384 (1992) 605. 

[83] S. Pirogov, Ya.G. Sinai, Theor. Math. Phys. Engl. Transl. 25 (1975) 1185; 26 
(1976) 39. 

[84] M. Zahradm'k, Commun. Math. Phys. 93 (1984) 559. 

[85] C. Borgs, J.Z. Imbrie, Commun. Math. Phys. 123 (1989) 305. 

[86] P.A. Slotte, P.C. Hemmer, J. Phys. C 17 (1984) 4645. 

[87] R.J. Baxter, J. Phys. C 6 (1973) L445. 

[88] T. Kihara, Y. Midzuno, T. Shizume, J. Phys. Soc. Jpn. 9 (1954) 681. 



28 



Table 1 

Exact results for the bulk energies Cq and e^, the latent heat Ae [87,88], and the 
values of V4{oo)\min and U2{oo)\max [16,27] for the q = 10, q = 8, q = 6, and q = 5 
Potts model in two dimensions. The last row refers to the R = 1 triangular SAP 
model. 





Model 


Co 


ed 


Ae 


V4{oo)\jnin 


U2{oo)\max 


q 


= 10 Potts 


-1.66425 


-0.96820 


0.696 


0.5589 


1.0751 


Q 


= 8 Potts 


-1.59673 


-1.11037 


0.486 


0.6207 


1.0333 


q 


= 6 Potts 


-1.50875 


-1.30775 


0.201 


0.6598 


1.0051 


q 


= 5 Potts 


-1.47367 


-1.42075 


0.053 


0.6662 


1.0003 




Tr. SAP 


-1.61784(728) 


-1.48852(565) 


0.129(12) 


0.66435(33) 


1.00174(25) 




Fig. 1. Energy distributions at T/j for the R = 1 triangular SAP model for selected 
lattice sizes. 
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Fig. 2. (a) Plot of the barrier height AF{L)L~'^~^^. The dotted hne is an extrap- 
olation to L — >^ oo, giving a nonzero value for the surface tension of the order of 
0.00539(12). The inset is used as a guide for the reader, (b) The energy minima 
eo(L) and e^iL) and their difference as a function of the inverse linear size. The 
dotted lines are linear fits indicating the values of the bulk energies Cq and e^, and 
that of the latent heat Ae at L = cxd. 




Fig. 3. Temperature-dependence of the ratios N{Vi,V2; P) (a) and Nx=2{Vi,V2; P) 
(b) (Eqs. (24) and (30), respectively). In panel (a) all pairs are shown, i.e. (Li, 
L2)=(40, 80), (50, 100), (60, 120), (70, 140), (80, 160), (100, 200), and (120, 240), 
while in panel (b) only the pairs (Li, -L2)=(80, 160), (100, 200), and (120, 240) are 
shown. 
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1.820 



Fig. 4. (a) Illustration of the crossing point corresponding to the peak of the num- 
ber-of-phases parameter. The behavior of the most probable energies of the pair 
of systems is also illustrated, (b) Behavior of the differences, (e) — e, of the two 
characteristic energies for each system of the chosen pair. For the larger lattice the 
graph corresponding to the double of the difference (e) — e is also shown in order to 
facilitate the illustration of the crossing relationship (31) for A = 2 (a = 4). 
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Fig. 5. (a) Simultaneous fittings for the shifting of all pseudocritical temperatures 
Tj* using an correction, (b) The same with panel (a) but with an correction. 
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Fig. 6. (a) Illustration of the behavior of the three peaks of the reduced num- 
ber-of-phases parameter defined in in Eq. (30) for A = 2. (b) Behavior of six finite- 
size transition points and comparison with our estimate (solid line) for the transition 
point and also with the estimate Tc = 1.8078 of Ref. [14] (dashed line). 




Fig. 7. FSS behavior of the specific heat peaks in two lattice ranges: L = 30 — 140 
and L = 100 — 480 (inset). An correction has been added to the divergence. 
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Fig. 8. (a) FSS behavior of the specific heat data at the five pseudocritical temper- 
atures T* using an correction, (b) The same with panel (a) but with an 
correction. Note the large difference in the values. 




Fig. 9. (a) FSS behavior of V/i{L)\min- Two fitting attempts are applied as shown in 
the figure, is much smaller for the case of an correction, (b) FSS of V4^{L)\min, 
now against L~^. The line corresponds to the fitting formulae shown in panel (b) 
and gives an estimate Vi{oo)\jnin = 0.66465(25). 
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Fig. 10. (a) The same as in Fig. 9 for U2{L)\max- (b) FSS of U2iL)\max, against L~^. 
The sohd Hne corresponds to the fitting formulae shown in panel (b) and produces 
the estimate U2{oo)\jnax = 1.00163(15). 




Fig. 11. (a) FSS behavior of the susceptibility data at the five pseudocritical temper- 
atures T*. The solid lines are simultaneous fitting attempts, using an correction, 
(b) The same with panel (a) but with a correction of the order of L~^. 
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